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Recent theoretical and imnierical evidence suggests that localization can survive the introduction 
of interactions in excited disordered many-body systems, giving rise to a so-called many-body lo- 
calization transition. This dynamical phase transition is relevant to questions of thermalization in 
extended quantum systems. It separates a many-body localized phase, in which localization prevents 
transport and thermalization, from a conducting ("ergodic") phase in which the usual assumptions of 
quantum statistical mechanics hold. Here, we present numerical evidence that many-body localiza- 
tion also occurs in models without disorder but rather a quasiperiodic potential. In one dimension, 
these systems already have a single-particle localization transition, and we show that this transition 
becomes a many-body localization transition upon the introduction of interactions. We also com- 
ment on possible relevance of our results to experimental studies of many-body dynamics of cold 
atoms and non-linear light in quasiperiodic potentials. 



I. INTRODUCTION 

In one-dimensional systems of non-interacting parti- 
cles, an arbitrarily weak disordered potential generically 
localizes all quantum eigenstatcs"'^'^. Such a system is 
always an insulator, with a vanishing conductivity in 
the thermodynamic limit. The question of how this pic- 
ture is modified by interactions remained unclear in the 
decades following Anderson's original work on localiza- 
tion^'^. Relatively recently, Basko, Aleiner, and Alt- 
shuler have argued that an interacting many-body system 
can undergo a so-called many-body localization (MBL) 
transition in the presence of quenched disorder. At low 
energy density and/or strong disorder, interactions are 
insufficient to thermalize the system, so the system re- 
mains a "perfect" insulator (i.e. with zero DC conduc- 
tivity despited being excited); at higher energy densi- 
ties and/or weaker disorder, the conductivity can become 
nonzero and the system thermalizes, leading to a con- 
ducting phase^'^. 

The MBL transition is rather unique for several rea- 
sons. First, in contrast to more conventional quantum 
phase transitions^, this is not a transition in the ground 
state. Instead, the MBL transition involves the local- 
ization of highly excited states of a many-body system, 
with finite energy density. This means that the transi- 
tion differs from most metal-insulator transitions, which 
are sharp only at zero temperature^. Furthermore, this 
MBL transition is of fundamental interest in the context 
of statistical mechanics. Local subsystems of interacting, 
many-body systems are generically expected to equili- 



brate with their surroundings, with statistical properties 
of these subsystems reaching thermal values after suffi- 
cient time. Studies of how this occurs in quantum sys- 
tems have led to the so-called eigenstate thermalization 
hypothesis (ETH), which states that individual eigen- 
states of the interacting quantum system already encode 
thermal distributions of local quantities''^". However, 
the many-body localized phase provides an example of a 
situation in which the ETH is false, and the ergodic hy- 
pothesis of quantum statistical mechanics is violated^^'^^. 
Since the work of Basko et al., these intriguing aspects of 
MBL have motivated many studies aimed at locating and 
understanding this transition in disordered systems^ 

On the other hand, it is important to note that single- 
particle localization does not require disorder. In 1980, 
Aubry and Andre studied a single-particle tight-binding 
model that omits disorder in favor of a potential that is 
periodic, but with a period that is incommensurate with 
the imderlying lattice^^. Much earlier. Harper considered 
this model as an effective one-dimensional description of 
an electron on a two-dimensional lattice with an irra- 
tional number of magnetic flux quanta per imit cell^"^. 
Aubry and Andre showed that such a quasiperiodic sys- 
tem exhibits a localization transition, even in one dimen- 
sion. The Aubry- Andre (AA) transition separates a weak 
potential phase, where all single-particle eigenstates are 
extended, from a strong potential phase where all eigen- 
states are localized. In the 1980s and 1990s, theorists 
continued to study this transition for its own peculiar- 
ities and because it mimics the situation in disordered 
systems in d > 3, where there is also a single-particle lo- 
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calization transition^*"^®. More recently, the AA model 
has been directly experimentally realized in cold atom 
experiments^^''^" and also in photonic waveguides'^^. The 
possibility of engineering quasiperiodic systems in the 
laboratory has inspired new theoretical and numerical 
work aimed at understanding the localization properties 
of such systems and how they differ from those with true 
disorder^^"^*. 



U A 



1. Statement of the problem and summary of the results 



In this paper, we ask whether there can be a MBL 
transition in an interacting extension of the AA model. 
More concretely, suppose we begin with a half-filled, one- 
dimensional system of fermions or hardcore bosons in a 
particular randomly chosen many-body Fock state, with 
some sites occupied and others empty. If the particles 
are allowed to hop and interact for a sufficiently long 
time, the standard expectation (e.g. in introductory 
treatments of statistical mechanics^^) is that the system 
should thermalize: that is, all microscopic states that are 
consistent with conservation laws should become equally 
likely and local observables should thereby assume some 
thermal distribution. Can this expectation be violated in 
the presence of a quasiperiodic potential, and if so, can 
this be traced to the persistence of localization even in 
the presence of interactions? 

The answer to both of these questions appears to be 
"yes." We use numerical simulations of unitary evolution 
of a many-body quasiperiodic system to measure three 
kinds of observables in the limit of very late times: the 
correlation between the initial and time-evolved parti- 
cle density profiles, the many-body participation ratio, 
and the Renyi entropy. Our observations are consistent 
with the existence of two regimes (phases) in the pa- 
rameter space of our model that differ qualitatively in 
ergodicity. At finite interparticle interaction strength 
u and large hopping 5, there exists a phase in which 
the usual assumptions of statistical mechanics appear to 
hold. The initial state evolves into a superposition of a 
finite fraction of the total number of possible configura- 
tions, and consequently, local observables approximately 
assume their thermal values. This is the many-body er- 
godic phase. However, at small hopping g, there is a 
phase in which particle transport away from the initial 
configuration is not strongly enhanced by interactions. 
The system explores only an exponentially small fraction 
of configuration space, and local observables do not even 
approximately thermalize. This is the many-body local- 
ized phase. Figure 1 presents a schematic illustration 
of the proposed phase diagram. Although interactions 
induce an expansion of the ergodic regime, the localized 
phase survives at finite w, and consequently, there is ev- 
idence for a quasiperiodic MBL transition. 
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FIG. 1: The proposed phase diagram for our interacting ex- 
tension of the Aubry-Andre model. Interactions convert the 
localized and extended phases of the AA model into many- 
body localized and ergodic phases and induce an expansion of 
the many-body ergodic phase. The phases of the interacting 
model differ qualitatively from their non-interacting counter- 
parts. The differences are explained in Section IV below. 



2. Organization of the paper 

We begin our study in Section II by introducing our 
interacting extension of the standard AA model. Since 
the MBL transition is a non-equilibrium phase transition, 
our goal is to understand the real-time dynamics of our 
many-body AA model. To simplify this task, we describe 
a method of modifying the dynamics of our model, such 
that numerical integration of the new dynamics is some- 
what easier than the original problem. In Section III, we 
introduce the quantities that we measure in our simula- 
tions and present the numerical results. Then, in Section 
IV, we argue that the data from Section III points to 
the existence of an MBL transition in our quasiperiodic 
system and extract information about the transition and 
phases. We conclude in Section V by summarizing our 
results, drawing connections to theory and experiment, 
and suggesting avenues for future extensions of our work. 

We relegate two exact diagonalization studies to Ap- 
pendix A. First, we examine the impact of our modified 
dynamics upon the single-particle and many-body prob- 
lems. Second, we study the many-body level statistics of 
the interacting model. We find evidence for a crossover 
between Poisson and Wigner-Dyson statistics, consistent 
with the usual expectation in the presence of a localiza- 
tion transition^ ^. 



II. MODEL AND METHODOLOGY 

In this section, we motivate and introduce our model 
and our numerical methodology for studying real-time 
dynamics. 
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A. The "Parent" Model 

We would like to consider one-dimensional lattice mod- 
els of the following general form: 

L-1 

^ =J2 ['^i^^ + (CjCj+i + c]+iCi) + Vnjfij+i (1) 

Here, cj is a fcrmion annihilation operator, and fij = CjCj 
is the corresponding fcrmion number operator. The three 
terms in the Hamiltonian (1) then correspond to an on- 
site potential, nearest-noiglibor hopping, and nearest- 
neighbor interaction respectively. For now, wc leave the 
boundary conditions unspecified. In ID, the Hamiltoni- 
ans (1) for hardcore bosons and fermions differ only in the 
matrix elements describing hopping over the boundary. 
With open boundary conditions, the Hamiltonians (and 
consequently all properties of the spectra) are identical. 

If wc set V = in the Hamiltonian (1) and take hj to 
be genuinely disordered, we recover the non-interacting 
Anderson Hamiltonian. If we then turn on a finite V ~ J, 
we obtain a model that is related to the spin models that 
have been studied in the context of MBL^^'^^. Alterna- 
tively, suppose we set F = again and take: 

hj=hcos{2TTkj + 5) (2) 

With a generic irrational wavenumber k and an arbitrary 

offset 5, we obtain the non- interacting AA modcl^^. For 
our purposes, we would like to use an incommensurate 
potential of the form (2), with h = l and 5 = ^ and u = 
^ left as tuning parameters to explore different phases 
of the model (1). 

Before proceeding, we should briefly review what is 
known about the single-particle AA model. This model 
is self-dual (for periodic boundary conditions and 5 = 0) 
and exactly solvable^^'^'"'^^. The self-duality can be re- 
alized by switching to Fourier space {Cj = e.^'^-'Cq) 
and then performing a rearrangement of the wavenum- 
bers q such that the real-space potential term looks like a 
nearest-neighbor hopping in Fourier space and vice versa. 
On a finite lattice of length L with periodic boundary 
conditions, such a rearrangement is possible whenever the 
wavenumber of the potential = j such that t and L are 
mutually prime. The duality construction reveals that, 
if the AA model has a transition, it must occur at g = ^ . 
In the thermodynamic limit, there is indeed a transition 
at this value for nearly all irrational wavenumbers k^^. 
When g > i, all single-particle eigenstates are spatially 
extended, and by duality, localized in momentum space; 
when g < i, all single-particle eigenstates are spatially 
localized, and by duality, extended in momentum space. 
Exactly at g = i, the eigenstates are multifractaP^. The 
spatially extended phase of the AA model is character- 
ized by ballistic, not diffusive, transport^^. Recently, Al- 
bert and Leboeuf have argued that localization in the 
AA model is a fundamentally more classical phenomenon 
than disorder-induced Anderson localization, and that 



the AA transition at g = ^ is most simply viewed as 
the classical trapping that occurs when the maximum 
eigenvalue of the kinetic (or hopping) term crosses the 
amplitude of the incommensurate potential^^. 

B. Numerical Methodology and Modification of 
the Quantum Dynamics 

Probing the MBL transition necessarily involves study- 
ing highly excited states of the system, and this precludes 
the application of much of the extensive machinery that 
has been developed for investigating low-energy physics. 
Consequently, several studies of MBL have resorted to 
exact diagonalization or other methods involving similar 
numerical cost^^'^^'^''. We too use a numerical method- 
ology that scales exponentially in the size of the system. 
However, in order to access longer evolution times in 
larger lattices, we introduce a modification of the quan- 
tum dynamics. This modification is inspired by a scheme 
used previously by two of us in a study of classical spin 
chains^^. There, at any given time, either the even spins 
in the chain were allowed to evolve under the influence 
of the odd spins or vice versa. This provided access to 
late times that would have been too difficult to access by 
direct integration of the standard classical equations of 
motion. 

By analogy, we propose allowing hopping on each bond 
in turn. At any given time, the instantaneous Hamilto- 
nian looks like: 

L-1 

Hm = LamJ(c\nCm+l + c\^+iCm) + ^ [hjUj + Vfljhj+i] 

J=0 

(3) 

We will specify the value of Om in Section II. C below, 
where we discuss our choice of boundary conditions. The 
state of the system is allowed to evolve under this Hamil- 
tonian for a time and this evolution can be imple- 
mented by applying the unitary operator: 

Um = exp (^-i^Hm^ (4) 

One full time-step of duration At consists of cycling 
through all the bonds: 

L-1 

U{At) =l[Um (5) 

m=0 

Note that, in (3), the hopping is enhanced by L because 
the hopping on any given bond is activated only once per 
cycle, while the potential and interaction terms always 
act. Therefore, the factor of L ensures that the average 
Hamiltonian over a time At has the form (1). The ad- 
vantage of employing the modified dynamics is that the 
Hm only couple pairs of configurations, so preparing the 
Um reduces to exponentiating order Vh two-by-two ma- 
trices, where Vh is the size of the Hilbert space. This is 
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a simpler task than exponentiating the original Hamilto- 
nian (1), which can take time cubic in Vh- This scheme 
only constitutes a polynomial speedup over exact diag- 
onalization, but that speedup can increase the range of 
accessible lattice sizes by a few sites. 

The modified dynamics raise several important is- 
sues that should be discussed^°. The periodic time- 
dependence of the Hamiltonian induces so-called "multi- 
photon" (or "energy umklapp") transitions between 
states of the "parent" model (1) that differ in energy by 
coh = reducing energy conservation to quasienergy 
conservation modulo ojh- Wc need to question whether 
this destroys the physics of interest: does the single- 
particle Aubry- Andre transition survive, or do the multi- 
photon processes destroy the localized phase? 

We take up this question in Appendix A, where we 
present a Floquet analysis of the singlci-particle and 
many-body problems. We find that, for sufficiently 
small At, the universal behavior of the single-particle AA 
model is presca'vcxi. At larger At, multi-photon processcis 
can strongly mix eigenstates of the single-particle parent 
model, increasing the single-particle density-of-states and 
destroying the AA transition. In the spirit of the earlier 
referenced work on classical spin chains^^, our perspec- 
tive in this paper is to identify whether MBL can occm- 
in a model qualitatively similar to our parent model (1). 
Therefore, to explore dynamics on long time scales, we 
avoid destroying the single-particle transition, but still 
choose At to be quite large within that constraint. 

In Appendix A, we also examine the consequences 
of our choice of At for the quasienergy spectrum of 
the many-body model. Our results suggest that multi- 
photon processes do not, in fact, strongly modify the par- 
ent model's spectrum for much of the parameter range 
that we explore in this paper^^. This means that partial 
energy conservation persists in our simulations despite 
the introduction of a time-dependent Hamiltonian, and 
we need to keep this fact in mind when we analyze our 
numerical data below. 

Finally, we note in passing that several recent stud- 
ies have focused on the localization properties of time- 
dependent models'*^^'*'^, including one on the quasiperi- 
odic Harper model"'^, but that the intricate details of this 
topic are somewhat peripheral to our main focus. 



C. Details of the Numerical Calculations 

In studies of the ID AA model, it is conventional to 
approach the thermodynamic limit by choosing lattice 
sizes according to the Fibonacci series (L = . . . 5, 8, 13, 
21, 34 . . .) and wavenumbers for the potential (2) as ratios 
of successive terms in the series^^. These values of k 
respect periodic boundary conditions while converging to 

the inverse of the golden ratio ^ = 0.618033 For any 

finite lattice, the potential is only commensurate with 
the entire lattice (since successive terms in the Fibonacci 
series are mutually prime), and the duality mapping of 
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TABLE I: For the various simulated lattice sizes L, the parti- 
cle number A'^, the configuration space size Vh, and the num- 
ber of samples used in the numerics. 

the AA model is always exactly preserved. 

For our purposes, this approach is unfortunately re- 
strictive, because lattice sizes above L = 21 arc inac- 
cessible. We have attempted to circumvent this limita- 
tion in several ways. For instance, we initially attempted 
to supplement the usual Fibonacci series with alternate 
Fibonacci-like series beginning with different seeds (e.g. 
L = ... 4, 7, 11, 18, 29 . . .). The ratio of successive terms 
of any such series converges to j-. Ultimately however, 
wc found empirically that finite-size effects are least prob- 
lematic if we use exclusively even L, always keep the 
wavenumber of the potential fixed at fc = ^, and set: 

Om = 1 - Srn,L-l (6) 

in equation (3), thereby forbidding hopping over the 
boundary^^. Note that, with these boundary conditions, 

our model describes hardcore bosons as well as fermions. 
The bosonic language maintains closer contact with cold 
atom experiments^^; the fermionic language is more in 
keeping with the MBL literature"'''^^. 

Using the approach described above, wc have simu- 
lated systems up to size L = 20 at half-filling. Our 
simulations always begin with a randomly chosen con- 
figuration (or Fock) state, so that the initial state has no 
entanglement across any spatial bond in the lattice, i.e. 
each particle is at some site with probability 1. Except 
in the exact diagonalization studies of Appendix A, we 
always set At = 1. Wc integrate out to tj = 9999 and 
ultimately average the evolution of measurable quantities 
over several samples, where a sample is specified by the 
choice of the initial configuration and offset phase to the 
potential (2). The sample counts used in the numerics 
are provided in Table I. 

III. NUMERICAL MEASUREMENTS 

We now introduce the quantities that we measure to 
characterize the different regimes of our model. We 
also present the numerical data along with some qual- 
itative remarks about the observed behavior. However, 
we largely defer quantitative phenomenology and mod- 
elling of the data to Section IV. 
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A. Temporal Autocorrelation Function 

One signature of localization is the system's retention 
of memory of its initial state. Note that we simulate 
the reversible evolution of a closed system. This means 
that the quantum state of the entire lattice retains full 
memory of its past, but we may still ask if the information 
needed to deduce the initial state is preserved locally or if 
it propagates to distant parts of the system. A diagnostic 
measure with which to pose this "local memory" question 
is the temporal autocorrelator of site j: 



(a) 



u = 0.32,9 = 0.1 



X,W^(2(n,)(t)-l)(2(n,)(0)-l) 



(7) 



Here, the angular brackets refer to an expectation value 
in the quantum state. This single-site autocorrelator may 
be averaged over sites and then over samples (as defined 
in Section II. C) to obtain: 



x{t;L) 



L-l 



(8) 



The sample average is indicated here with the large 
square brackets. Typically, to reduce the effects of noise, 
we also average over a few time steps within each sam- 
ple (i.e. perform time binning) before taking the sample 
average. 

We can discriminate three qualitatively different be- 
haviors of X '^s. t in our interacting model. Figure 2 
shows examples of each of these behaviors at interac- 
tion strength u — 0.32. Panel (a) is characteristic of 
the low g regime, where the autocorrelator stays invari- 
ant over several orders-of-magnitude of time, and there is 
no statistically significant difference between time series 
for different L. At higher g, as in panel (b), the time 
series show approximately power-law decay culminating 
in saturation to a late-time asymptote. For the largest 
systems, the power law is roughly consistent with the dif- 
fusive expectation of t~ 5 decay. The late-time asymptote 
decays with L (as expected from energy conservation^^ ) 
suggesting that the power-law decay may continue indef- 
initely in the thermodynamic limit. Surprisingly, at still 
larger there is a third behavior, exemplified by panel 
(c). For the largest lattice sizes, the power-law era is not 
followed by saturation but by an extremely rapid decay. 
The rapid decay is most evident in the large g, large u 
regime, where the energy density of the parent model (1) 
is relatively large. This implies that this behavior might 
be tied to the multi-photon processes induced by peri- 
odic modulation of the Hamiltonian; correspondingly, it 
also implies that, for fixed g and u, we might be able to 
induce the appearance of the rapid decay by increasing 
At. We have tested this numerically, and the results sup- 
port the connection to the energy non-conserving multi- 
photon processes. This suggests that there are only two 
distinct regimes of the parent model represented in Fig- 
ure 2, differentiated by the L dependence of the asymp- 
totic value of the autocorrelator. We will proceed under 
this working assumption. 




(b) 



u = 0.32, g = 0.5 



(c) 




FIG. 2: Three characteristic time series for the temporal auto- 
correlator with u = 0.32 and At = 1. In each panel, we show 
time series for a particular value of the hopping g. Only a 
few representative error bars are displayed in each time series. 
The legend refers to dilTerent lattice sizes L. The reference 
lines in panels (b) and (c) show diffusive t~^ decay. 



The difference between these two regimes is brought 
out more clearly in Figure 3. We focus on a late time t = 
itost and probe x(^tcst; L) as a function of g for different 
lattice sizes. Panels (a)-(c) show data for u = 0, 0.04, 
and 0.64 respectively. All the panels show a "splaying" 
point of the x vs. L curves, separating a high g regime 
in which x(^tost; L) decays with L from a low g regime in 
which it does not. The value of g at this feature decreases 
monotonically with u. Most importantly, in each case, 
this value is robust to changing itost ; if we halve ttest from 
the value that appears in Figure 3, the feature appears 
at approximately the same value of g. This property 
of the data is very fortunate: in Section IV. C below, 
we will use the splaying feature in these plots to put a 
numerical lower bound on the transition value of g for 
different interaction strengths. Since time scales get very 
long near the transition, it is difficult to simulate out to 
convergence in this regime. Nevertheless, the fact that 
the value of g at the splaying feature remains fixed in time 
implies that we can deduce the phase structure from our 
finite-time observations. 
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(a) 



u = 0, t. =9980-9999 



(b) 



(c) 




FIG. 3: The value of x in the latest time bin {t = 
9980 ... 9999) plotted against g. In panels (a)-(c), u = Q, 
0.04, and 0.64 respectively. The legend refers to different lat- 
tice sizes L. 



B. Normalized participation ratio 

One of the commonly used diagnostics for studying 
single-particle localization is the inverse participation ra- 
tio (IPR). This quantity is intended to probe whether 
quantum states explore the entire volume of the system 
and is often defined as the sum over sites of the ampli- 
tude of the state to the fourth power: Y^- IV'il'*- Typi- 
cally, the IPR is independent of system size in a localized 
phase and decays as the inverse of system volume in an 
extended phase. We now describe how this quantity can 
be fruitfully exploited in the many-body context. Let 
c denote some specific configuration of N particles in L 
sites. Then, we can write the state of the system in the 
configuration basis as: 



\^{t))=Y.i^c{t)\c) 

The configuration-basis IPR is simply: 



P{t-L) = 



(9) 



(10) 



where the square brackets, as usual, denote a sample av- 
erage. 

Identifying a many-body localized phase using the IPR 
(10) is more tricky than in the single-particle context. 
This is because a localized many-body system may still 
occupy an exponentially large (in L) number of configu- 
rations. A cartoon picture of this emerges if we consider 
a system of distinguishable particles evolving from an ini- 
tial configuration state in the non-interacting Anderson 
Hamlitonian. The number of configurations explored by 
the system scales as where ^ is the localization length, 
and the IPR scales as the inverse of this number: . 
At half-filling N = so the IPR decreases exponen- 
tially with L. Meanwhile, in a perfectly ergodic phase, 
the system explores all possible configurations that are 
consistent with particle number conservation. The num- 
ber of such configurations is also exponential in system 
size. However, a sharp distinction between ergodic and 
localized behavior appears in the normalized participa- 
tion ratio (NPR): 



1 



P{t\L)VH 



(11) 



We expect rj{t; L) to be independent of L at late times 
in the ergodic phase. In a many-body localized phase, 
although Pit; L) would decay exponentially, it would do 
so more slowly than and rj{t\L) would decay expo- 
nentially in L as well. 

In Figure 4, we plot ri{ttcst \ L) vs. g for it ^ 0, 0.04, and 
0.64. The figure reveals an important difference between 
the non-interacting and interacting models. At low g, 
both with and without interactions, rj decays exponen- 
tially with L: 



rj (X exp(— kX) 



(12) 



with K > 0. More surprisingly, 77 also decays with L at 
large g in the non-interacting case; all that happens is 
that K becomes essentially independent of g. With even 
small interactions however, rj becomes system-size inde- 
pendent in the large g regime, following our ansatz for an 
ergodic phase. We bring out this point more clearly in 
Figure 5, in which we extract estimates for the decay co- 
efficient K for various values of the interaction strength. 
Thus, the extended phase of the non-interacting AA 
model appears to be a special, non-ergodic limit. 

Before proceeding, we should caution that, in panels 
(b) and (c) of Figure 4, the collapse at high g looks very 
appealing because of the use of a semilog plot and would 
not be so striking on a normal scale. The axes have 
been chosen to highlight the exponential scaling at low 
g, which would not be as apparent if we simply plotted 
Tj vs. g. However, regarding the absence of perfect col- 
lapse at high g, note that the raw data for the IPR differ 
by several orders-of-magnitude for different values of the 
lattice size L. Given this, the coincidence of the order-of- 
magnitude of for different values of L is already a good 
indication of the proposed scaling, and some corrections 
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(b) 



u = 0.04, t. . = 9980-9999 




(c) 
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FIG. 4: The value of r\ in the latest time bin it = 
9980 ... 9999) plotted against g. In panels (a)-(c), u = 0, 
0.04, and 0.64 respectively. The legend refers to different lat- 
tice sizes L. See equation (11) for the definition of r\. In the 
ergodic phase -r] « 0.5. 



t. =9980-9999 




FIG. 5: Estimates of k from a fit of oc e in the latest 
time bin (tbin = 9980 - 9999). The legend refers to diflterent 
values of the interaction strength u. 



to this scaling should be expected given the modest ac- 
cessible system sizes. 



C. Renyi Entanglement Entropy 

Unlike the normalized participation ratio, which pro- 
vides a global characterization of the time-evolved state, 
bipartite entanglement is arguably a better proxy for 
whether a part of the system can act as a good heat 
bath for the rest. In the many-body ergodic phase, we 
expect the bipartite entanglement entropy to be a faith- 
ful reflection of the thermodynamic entropy. This im- 
plies an extensive entropy, pinned to its thermal infinite 
temperature value throughout the phase^^ In contrast, 
in the many-body localized phase, we expect an exten- 
sive but subthermal entanglement entropy. This expec- 
tation is consistent with the results of three recent papers 
that focus on the behavior of entanglement measures in 
the many-body localized phase of the disordered prob- 
\%va^^'^'^''^^ . These papers also study the time dependence 
of the entropy beginning from an unentangled product 
state. In the many-body localized phase, this growth is 
found to be slow, generically logarithmic in time. Since 
our model lacks disorder altogether, it may be interest- 
ing to explore the entanglement dynamics here as well. 
In what follows, we comment on the dynamics, but we 
primarily use the late-time entanglement entropy as yet 
another tool to help distinguish between the many-body 
locahzed and ergodic phases. 

Let subsystem A refer to lattice sites 0, 1, ... ^ — 1, 
and let subsystem B refer to the remaining sites in the 
chain. We can compute the reduced density matrix of 
subsystem A by beginning with the full density matrix 
/5(t) = l^'(t)) (^'(01 ^iid tracing out the degrees of free- 
dom associated with subsystem B: 



(13) 



The sample-averaged order-2 Renyi entropy of subsystem 
A is then given by: 



52(i;L)^ [-l0g2 {TTA{pAitf})\ 



(14) 



Both 5*2 and the standard von Neumann entropy are ex- 
pected to attain the same values in the ergodic phase; 
we choose to focus on the former to save on the compu- 
tational cost of diagonalizing the reduced density matrix 

(13). ^ ^ ' 

Our first task is to examine whether the putative local- 
ized phase of our model exhibits the same behavior that 
was observed with tDMRG^"^'^^. In panel (a) of Figure 
6, we focus on a low value of g and plot ^2 vs. ln(i) 
for i = 10 lattices. At very early times, the time series 
all tend to coincide, reflecting the formation of short- 
range entanglement at the cut between the subsystems. 
Afterwards, the non-interacting time series saturates for 
several orders-of-magnitude of time, while the interacting 
time series show behavior that is consistent with logarith- 
mic growth. In order to clearly establish the saturation 
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ln(t) 

FIG. 6: Example time series of the Renyi entropy for two val- 
ues of the tuning parameter g. The legend refers to different 
values of the interaction strength u. Panel (a) shows data for 
L = 10 lattices at p = 0.2. Panel (b) shows data for L = 20 
lattices at p = 1.1. In the localized regime, we need to use 
smaller lattices to see convergence Renyi entropy. 

that follows the slow growth, we have had to focus on 
small lattices. Panel (b) of Figure 6 shows data for large 
g. Here, the most striking difference between the non- 
interacting and interacting models lies in the saturation 
value of the entropy: the interacting model is substan- 
tially more entangled, but the saturation value does not 
appear to depend on the value of u. We will see below 
that this is another indication that thermalization only 
occurs in the interacting, large g regime. 

Figure 7 shows late-time values of the Renyi entropy 
density plotted against the tuning parameter g. We first 
focus on the high g regime. In panel (a), u = 0, and 
S2{ttcsuL) cx L for large g. However, the entropy den- 
sity is less than ^, which is the thermal result when the 
system has ergodic access to all configurations consistent 
with particle number conservation. The situation is dra- 
matically different in panels (b) and (c), where u — 0.04 
and 0.64 respectively. At high g, the entropy actually 
looks superextensive. This is just a finite-size effect, be- 
cause the entropy is well fit to a linear growth of the 
form: 

5'2(<tcst; L) = mL - Sdci (15) 

where S'dcf is a constant deficit, typically around 1.15 — 
1.3. In Figure 8, we show that the slope m « ^ at large g 
in the interacting problem. This implies that the entropy 
is thermal in the L — >■ oo limit, where the deficit S'dcf is 
negligible. 

Now, we turn to the low g regime, where all entangle- 
ment formation is local to the cut between the subsys- 



tems. Without interactions, the off-diagonal elements in 
the reduced density matrix (13) typically contain only a 
few frequencies originating from localized single particle 
orbitals immediately adjacent to the cut. The number 
of relevant orbitals is finite in L. As a result, the off- 
diagonal elements cannot fully vanish, and the reduced 
density matrix never thermalizes. The resulting entan- 
glement entropy is independent of L as shown in the inset 
of panel (a). In the interacting problem, while the or- 
bitals immediately adjacent to the cut still have roughly 
the same frequencies, the "spectral drift" (i.e. the spread 
of these lines due to sensitivity to the configuration of dis- 
tant particles) allows for a much larger number of distinct 
and mutually incoherent contributions to offdiagonal ele- 
ments of the reduced density matrix. These off-diagonal 
elements can dephase more efficiently, leading to a par- 
tial thermalization. This is the mechanism that likely 
underlies the extensive but subthermal entropy observed 
by Bardarson et al.^^. For small L, our numerical results 
in the low g regime agree well with this expectation. For 
larger L, the slow dynamics of the entropy formation 
makes it difficult to observe saturation, both in our work 
and in the tDMRG of Bardarson et al. 

If the entropy eventually becomes extensive for all L, 
then the "crossing" feature that is present in panels (b) 
and (c) of Figure 7 would become a "splaying" feature, 
with the entropy density independent of L at small g. In 
any case, an interesting property of the data is that the 
values of g at the crossing features of the S'2(itcst; L) vs. 
g plots are consistent with the locations of the splaying 
features in the corresponding xiUcst', L) vs. g plots of 
Figure 3. This seems to be the case for all u. Thus, these 
features may be useful in locating the transition. 



IV. DISCUSSION 

Above, we presented numerical evidence that our in- 
teracting AA model contains two regimes that show qual- 
itatively distinct behavior of the autocorrelator, normal- 
ized participation ratio, and Renyi entropy. Next, we 
will propose and characterize model quantum states that 
qualitatively (and sometimes quantitatively) reproduce 
the numerically observed late-time behavior in the two 
regimes. These model states expose more clearly why 
the two regimes of our model are appropriately identified 
as many-body ergodic and localized phases. To conclude 
this section, we will also use our numerical data to trace 
the phase boundary separating these phases, thereby mo- 
tivating the phase diagram in Figure 1. 



A. The Many-Body Ergodic Phase 

To model the behavior of the putative ergodic phase, 
we begin by writing down a generic model state in the 
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(a) 



u = 0, t = 9999 



(b) 



(c) 




0.2 0.4 0.6 0.8 1 1.2 1.4 



FIG. 7: The value of ^ at t = 9999 plotted against g. In 
panels (a)-(c), m = 0, 0.04, and 0.64 respectively. The legend 
refers to different lattice sizes L. In panel (a), the inset plot 
shows 52 vs. g in the low g regime. In panels (b) and (c), the 
for low L in the low g regime. 



insets show ^ vs 



t = 9999 




FIG. 8: The estimated slope of 5*2 vs. L at late times as 
a function of g. The legend refers to different values of the 
interaction strength u. 



configuration basis: 
1$) 



{c} 



E E 

n=0 {ca,cb} 



An) 



Sn) (n) 



(16) 



Here, the c refer to configurations of the full chain, 
whereas the ca and cb refer to configurations of the sub- 
systems A and B, as defined in Section III.C above. The 
superscripts on the configurations and expansion coeffi- 
cients refer to the number of particles in subsystem A. 
Writing the state in terms of the subsystem configura- 
tions will be useful shortly, but for now we focus on the 
statistical properties of the amplitude (pc We assume 
that this amplitude is distributed as a complex Gaussian 
random variable: 



1 



27rcr' 



■ exp 



2(t2 



(17) 



Within this distribution, (l^p) = 2ct2 and = Sct'*. 

From these average values, it is possible to deduce that: 

^ (18) 



for normalization and that the IPR is P$ = This, in 



turn, implies: 



(19) 



This result is reproduced quantitatively in the numerics 
in Figure 4. 

Next, suppose we compute the reduced density matrix 
of subsystem A in the state |$): 



PA 



E E 

n {ca,Cj^,,cb} 



/,*(") J,(«) 
^AB Va'B 



(20) 



To find the Renyi entropy, we need to compute the trace 
of the square of this operator: 



TrA{p\} 



E E 



f^AB 'l'A'B't>AB'(t^A'B' 



•} 



When we average over our distribution of amplitudes 
(17), only the coherent terms survive '''2: 



TrA{p\} 



E E 

+E E 

n {ca,ca'-^ 

-E E 

n {ca,cb} 



{\&\^^a^b\') 



(22) 



The final term accounts for the double counting of terms 
where ca = ca' and cb = cb' simultaneously. We now 
introduce the notation: 



p\ 

Q\{P~Q)\ 



(23) 



10 



and evaluate the expectation values in equation (21) to 
obtain: 



of the form: 



(24) 



Finally, using a Stirling approximation to the combina- 
tion function and a saddle-point approximation for the 
sum, we find the entropy: 



52,*«§-l0g,(-| 



L 



1.2 



(25) 



This is the same form observed in the numerics (15), and 
the deficit 5def lies in the observed range. Asymptotically 
in L, the entropy (25) is maximal^ and this is exactly the 
expected behavior when the particle number thcrmalizes. 

There is an important caveat to note here: we have 
argued above that, if multi-photon processes do not com- 
pletely destroy energy conservation, then this can lead to 
relic autocorrelations at late times. This implies that the 
assumption of independent random amplitudes cannot be 
exactly correct on a finite lattice. However, the relic au- 
tocorrelations decay with L, suggesting that our assump- 
tions get better as the system size grows. Therefore, in 
the thermodynamic limit, this phase is truly thermal. 



B. The Many-Body Localized Phcise 

Our model for the time-evolved state in the localized 
regime is founded upon the following intuition: there ex- 
ists a length scale ^, which is analogous to the single- 
particle localization length and beyond which particles 
arc unlikely to stray from their positions in the initial 
state. Then, if we partition our lattice of length L into 
blocks of size ^, exchange of particles between blocks 
is less important than rearrangements of the particles 
within each block. Consequently, the total number of 
configurations accessed by the state of the full system 
is approximately the product of the number of configu- 
rations accessed within each block. This multiplicative 
assumption should be very safe in a localized phase. We 
additionally assume that, within each block, the dynam- 
ics completely scramble the particle configuration. If a 
certain block of length ^ contains n particles in the initial 
state, then the time-evolved state contains equal ampli- 
tude for each of the possible ways of arranging n particles 
in those ^ sites. In keeping with our numerical protocol, 
we randomly select the initial state from the space of all 
possible Fock states of a certain global particle number. 
Then, a block of ^ sites contains n particles with proba- 
bility: 



w{£,,n) = 



2« 



l + O 



2 \ 1 



(26) 



We will consider the limit L » ^ 1, where we can 

approximate the probability by the first term. The as- 
sumptions proposed above motivate writing down a state 



|A) 



1 



{ci,...cl} 



(^ci,...CL^ ci,...C|^ (27) 



where the tilde on the sum indicates that it should only 

run over configurations that are consistent with the initial 
distribution of particles among the blocks. The factors z 
are complex phases which depend upon the configuration, 
and M is a normalization which is equal to the total 
number of configurations represented in the state |A). 

Before beginning our analysis of the state |A), we 
should note that, in contrast to our calculations in the 
ergodic phase, our goal in the localized regime will be 
to qualitatively tie the numerically observed large L be- 
havior to the existence of the length scale ^. Unfortu- 
nately, we cannot achieve the quantitative accuracy of 
the ergodic model state |$) with the localized toy- model 
described above. 

We begin by estimating the autocorrelator between the 
initial state and the model time-evolved state |A). A non- 
zero autocorrelator emerges, because each block is only 
at half-filling on average. Fluctuations away from half- 
filling (in either direction) yield a positive typical value of 
the autocorrelator within a block. Indicating an average 
over the distribution (26) with an overline, wc find the 
block value Xbiock ~ j • This is also the average value for 
the whole system when i » ^: 



XA 



(28) 



Next, to estimate the IPR, we need to compute the nor- 
malization factor M . We begin by estimating the number 
of explored configurations in each block. The average of 
the logarithm of the number of explored configurations 
within a block is: 

HMuocV) « In - ^ (29) 

Then, using InM « |^lnMbiock, we can estimate M itself 
as: 



(30) 

Using this normalization, we can estimate the NPR ?7a: 

1 



InryA 



2 1^^+^ 1^(1 



(31) 



This qualitatively agrees with the numerically observed 
behavior (12) up to subleading corrections, and in the 
large- L limit: 



2^ 



(32) 



Note that equations (28) and (32) imply a relationship 
between the scaling behaviors of x and k in the localized 
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regime. This relationship is not reflected in our numerical 
data, in part because we cannot truly attain the limit 
i » ^ 1. The numerically computed value of k, for 
example, can contain finite-size corrections of order 

or ^. Also, we must keep in mind that the state |A) 
is just a toy model that does not capture fine details of 
the time-evolved states in this regime. Thus, we must be 
content with reproducing the qualitative behavior of each 
measurable quantity individually, without expecting the 
relationships between these quantities in | A) to be exactly 
reproduced in the data. 

We now turn to the Renyi entropy, the quantity which 
most strikingly distinguishes between the non-interacting 
and interacting localized phases. To examine this quan- 
tity, we revert to partitioning the system in half, instead 
of into blocks of size ^. As long as ^ -C ^ , the assump- 
tions that we made above about the blocks of size ^ hold 
even better for the subsystems A and B. For example, we 
can assume that there are "explored sets" of Ma config- 
urations in subsystem A and Mb configurations in sub- 
system B respectively, with M = MaMb- We consider 
computing the reduced density matrix pA, exactly as in 
equation (20) above. If the off-diagonal elements of this 
density matrix remain perfectly phase-coherent, it can 
easily be shown that 5*2°^' = 0. In reality, there will be a 
local contribution to the entropj^ from particles straying 
over the cut between subsystems A and B. This mim- 
ics the situation in non-interacting localized phases. Al- 
ternatively, suppose that dephasing is sufficiently strong 
that we can proceed by analogy with the ergodic phase, 
beginning with equation (21) and keeping only coherent 
terms as in equation (22). Thereafter, the calculation for 
the model localized state |A) differs from the calculation 
for |<I>). We need to consider the statistics of the con- 
figm'ation probabilities |A^bP- For IA^^bP 7^ 0, we need 
the configurations on both subsystems to lie within their 
respective explored sets; this occurs in subsystem A, for 
example, with probability ■ This reasoning leads 

to the "dephased" entropy: 



ridp 



■l0g2 

■l0g2 
1 - 





1 




Ms ' 


^^- 




V\/m 


M J 











1 



MaM. 



(33) 



where we have additionally made the approximation that 
typically Ma ~ Mb « \/M. With only partial loss of co- 
herence, the entropy would lie between these two limiting 
cases: S^jV < a < 'S'2 a- Thus, dephasing alone, with- 
out additional particle transport, can induce an extensive 
entropy. 

Indeed, our numerics support the view that the main 
difference between the non-interacting and many-body 

localized phases is the amount of dephasing. There does 
not seem to be a qualitative difference in particle trans- 



u 


X 


At 


m 


0.04 


0.35 


0.45 


0.45 


0.16 


0.30 


0.40 


0.40 


0.32 


0.25 


0.40 


0.40 


0.64 


0.25 


0.40 


0.35 



TABLE II: Bounds or estimates of the transition value of gc 
at various values of u and based on various measured quanti- 
ties. See the text for the reasoning behind the estimates. All 
values carry implicit error bars of ±0.05 as that is the dis- 
cretization of our simulated values of g. This error bar should 
be interpreted, for instance, as the error on our estimate of 
the location of the maximum value of m. The error on our 
estimate of gc is, of course, much larger. 



port. The particle configuration stays trapped near its 
initial state, even with interactions, and the system does 
not thermalize. 



The Quasiperiodic Many-Body Localization 
Transition 



Estimating the location of the MBL transition is ex- 
tremely challenging. Given the numerically accessible 
lattice sizes, satisfying finite-size scaling analyses are dif- 
ficult to perform. Nevertheless, rough estimates have 
been made in the disordered problcm-'^'^''^^'"'^^'^"'^, so wc will 
now attempt to extract an approximate phase boundary 
for our model. 

We first consider the autocorrelator. Above, wc noted 
the "splaying" feature in the late-time plots of the auto- 
correlator vs. g. The value of g at this feature can be 
taken as a lower bound for the transition. For g slightly 
greater than this value, it is possible that x only decays 
with L because ^ > L for accessible lattice sizes. Consid- 
ering two lattice sizes {L = 16 and 20) and finding when 
their values of x deviate, we find the values reported in 
the first column of Table II. 

Next, we consider the fitting parameter k in equation 
(12). In Figure 5, we see that there is a region where 
K < for finite interaction strength. Since rj < I, finite- 
size effects are clearly dominating the estimate in this 
region. We can use the value of g where k is minimal to 
track how this region moves as u is varied. This yields 
the second column of the table. 

Finally, a similar approach can be applied to extract 
estimates of gc from the fits (15). There exists a region 
where m > 5, but this is mathematically inconsistent in 
the thermodynamic limit. Therefore, if we find the value 
of g that maximizes m, we can again estimate the location 
of the region dominated by finite-size effects, yielding the 
final column of Table II. 

The estimates of the transition value gc in Table II 
were obtained using data for the latest time that we sim- 
ulated (the time bin fbin = 9980 . . . 9999 for x and k and 
t = 9999 for m). However, we have also estimated gc for 
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data obtained at a half and a quarter of this integration 
time, finding consistent results. Thus, the general phase 
structure of the model is invariant to changing the obser- 
vation time, even though not all measurable quantities 
have converged to their asymptotic values. Consolidat- 
ing the information from the estimates in Table II, we 
propose that the phase diagram qualitatively resembles 
Figure 1. 

Before proceeding, it is worth noting that our rough 
estimates of the phase boundary do not make assump- 
tions regarding the character of the MBL transition (i.e. 
whether it is continuous or first order). In fact, some of 
our plots (e.g. panel (c) of Figure 7) hint at the possi- 
bility of a discontinuous change in 52 as a function of g 
in the thermodynamic limit. Wc are not aware of any 
results that rule out a first-order MBL transition, so we 
must keep this possibility in mind. 



V. CONCLUSION 

Recently, evidence has accumulated that Ander- 
son localization can survive the introduction of suffi- 
ciently weak interparticle interactions, giving rise to 
a many-body localization transition in disordered sys- 
tems"'''^'^^'^^'^^. The MBL transition appears to be a 
thermalization transition: in the proposed many-body 
localized phase, the fundamental assumption of statisti- 
cal mechanics breaks down, and the system fails to serve 
as its own heat bath^^'^^. We have presented numeri- 
cal evidence that this type of transition can also occur 
in systems lacking true disorder if they instead exhibit 
pseudodisorder in the form of a quasiperiodic potential. 

From one perspective, this may be an unsurprising 
claim. For g < 5 the localized single-particle eigenstates 
of the quasiperiodic Aubry- Andre model have the same 
qualitative structure as those of the Anderson model, so 
the effects of introducing interactions ought to be similar. 
By this reasoning, perhaps it is even possible to guess the 
phase structure of an interacting AA model using knowl- 
edge of an interacting Anderson model: we simply match 
lines of the two phase diagrams that correspond to the 
same non-interacting, single-particle localization length. 

However, this perspective misses important effects in 
all regions of the phase diagram. Most obviously, the AA 
model has a transition at w = 0, and it is interesting to 
see how this transition gets modified as it presumably 
evolves into the MBL transition at finite u. It is also 
important to remember that quasiperiodic potentials are 
completely spatially correlated. This means that the AA 
model lacks rare-regions (Griffiths) effects, and this may 
have subtle consequences for the dynamics. Finally, the 
AA model contains a phase that is absent in the one- 
dimensional Anderson model, the .9 > 5 extended phase, 
and we have seen above that interactions have a profound 
effect upon this regime. 

Understanding MBL in the quasiperiodic context is es- 
pecially pertinent given the current experimental situa- 



tion. Some experiments that probe localization physics 
in cold atom systems use quasiperiodic potentials, con- 
structed from the superposition of incommensurate opti- 
cal lattices, in place of genuine disorder. The group of In- 
guscio, in particular, has recently explored particle trans- 
port for interacting bosons within this setup^^'^°. Mean- 
while, the A A model has also been realized in photonic 
waveguides, and experimentalists have studied the effects 
of weak interactions on light propagation through these 
systems. They have also investigated "quantum walks" 
of two interacting photons in disordered waveguides^^'^^. 
This protocol resembles the one we have implemented nu- 
merically, so similar physics may arise. Finally, we note 
that Basko et al. have predicted experimental manifes- 
tations of MBL in solid-state materials. In such systems, 
there is always coupling to a phononic bath, so the MBL 
transition is expected to become a crossover that nev- 
ertheless retains interesting manifestations of the MBL 
phenomena*^. Whether there exist quasiperiodic solid- 
state systems to which the predictions of Basko et al. 
apply remains to be understood. 

Given the current experimental relevance of localiza- 
tion phenomena in quasiperiodic systems, we hope that 
our study will motivate further attempts to understand 
these issues. Unfortunately, our ability to definitively 
identify and analyze the MBL transition is limited by 
the modest lattice sizes and evolution times that we can 
simulate. Vosk and Altman recently developed a strong- 
disorder renormalization group for dynamics in the disor- 
dered problem^", but the reliability of such an approach 
in the quasiperiodic c:ontext is unclear. A time-dependent 
density matrix renormalization (tDMRG) group study 
of this problem would be a valuable next step. Tezuka 
and Garcfa-Garcfa have published tDMRG results on lo- 
calization in an interacting AA model, but their focus 
was not on the thermalization questions of many-body 
localization^^. It would be worthwhile to pose these ques- 
tions using a methodology that allows access to much 
larger lattices. However, even tDMRG may have diffi- 
culty capturing the highly-entangled ergodic phase^^'^^, 
so an effective numerical approach for definitively char- 
acterizing the transition remains elusive. 
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Appendix A: Exact Diagonalization Results for the 
Single-Particle and Many-Body Problems 

This appendix collects exact diagonalization results 
that supplement the real-time dynamics study in the 
main body of the paper. 



1. Floquet Analysis of the Modified Dynamics 

The goal of the first part of this appendix is to exam- 
ine the consequences of the modifications to the quantum 
dynamics described in Section II. B above. We first ver- 
ify that the AA transition survives by diagonalizing the 
single-particle AA Hamiltonian (i.e. the Hamiltonian (1) 
with u = = 0) and the single-particle unitary evolu- 
tion operators (5) for various choices of the time step At. 
Subsequently, we employ the same approach to examine 
how varying At impacts the quasienergy spectrum of the 
interacting, many-body model. 

a. Robustness of the Smgle-Particle Aubry-Andre 
Transition 

To study the single-particle transition, we focus on the 
inverse participation ratio: 

PsAa-.L)-^ (^El^^l'j (Ai) 

Here, ipj denotes the amplitude of the wave function at 
site j of an L site lattice. We enclose the sum in equation 
(Al) in parentheses to indicate important differences in 
the averaging procedure with respect to the many-body 
inverse participation ratio (10). In the many-body case, 
we computed the IPR as a sum over configurations in the 
quantum state at a particular time in the real-time evo- 
lution. Then, we averaged over samples, where a sample 
was specified by a choice of the offset phase to the po- 
tential (2) and an initial configuration. Throughout this 
appendix, we instead specify a "sample" solely by the off- 
set phase S, and we average over eigenstates within each 
sample before averaging over samples. 

As noted previously, the usual AA model has a tran- 
sition that must occur, by duality, a,t gc = \- Near the 
transition, the localization length is known to diverge 
with exponent v — 1^^. Our exact diagonalization re- 
sults indicate that, at the transition, Psp{gc, L) ^ L^^. 




(g-0.45)L 

FIG. 9: Collapse of single-particle IPR vs. g, using the scaling 
hypothesis (A2). The legend refers to different lattice sizes L. 
In panel (a), we show data for the usual AA Hamiltonian (1). 
In panel (b), we show data obtained from diagonalizing the 
unitary evolution operator for one time step in the modified 
dynamics (5). We use potential wavenumber fc = ^ and 50 
samples for all lattice sizes. 

Hence, we can make the following scaling hypothesis for 
the IPR: 

Psp = L-^-f{{g-g,)L) (A2) 

In panel (a) of Figure 9, we show that we can use this 
scaling hypothesis to collapse data for the standard AA 
model. We show data for i = 8 to L = 128, with poten- 
tial wavenumber fc = ^ and open boundary conditions. 
For all lattice sizes, we average over 50 samples. 

To establish the stability of the AA transition to the 
modified dynamics, we must ask: can the IPR obtained 
from diagonalizing the unitary evolution operators (5) be 
described using the scaling hypothesis (A2)? Panel (b) 
of Figure 9 shows that this is indeed the case for At = 1. 
The only parameter that needs to be changed is gc, which 
decreases slightly as At is raised. This implies that there 
is a transition in the Floquet spectrum of the system that 
can be tuned by varying At. It would be a worthwhile 
exercise to map out the phase diagram of this single- 
particle problem in the {g, At) plane. We leave this for 
future work. 



b. Properties of ttie Many-Body Quasienergy Spectrum 

We now turn our attention back to the effects of the 
modified dynamics upon the full, many-body model. 
In Section II. B above, we emphasized that our time- 
dependent model lacks energy conservation, with multi- 
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FIG. 10: The density-of-states vs. quasienergy for L = 12 
systems at half-filling with interaction strength u = 0.16. The 
legend refers to different values of At; the time- independent, 
parent model is referred to as "PM." In panels (a)-(c), g = 
0.25, 0.4, and 0.9 respectively. 



photon processes inducing transitions between states of 
the "parent" model (1) that differ in energy by ujh — 
In this part of the appendix, we will examine how vary- 
ing At impacts the quasienergy spectrum of the time- 
dependent model, using the approach that we applied to 
the single-particle case above: we diagonalize the time- 
independent Hamiltonian as well as the unitary evolution 
operator for one time step of the time-dependent model. 

In Figure 10, we plot the density-of-states d{At, E) 
in quasienergy space of the parent model and time- 
dependent models for different values of Ai. We focus 
on L = 12 systems at half-filling with fermions (or, since 
we continue to use the boundary conditions described 
in Section II. C, hardcore bosons). We fix the interac- 
tion strength to m = 0.16 and tune g to explore different 
regimes of the model. In panels (a)-(c), we plot data for 
g = 0.25, 0.4, and 0.9. According to Table II, these val- 
ues of g put the system in the localized phase, near the 
transition, and in the ergodic phase respectively. 

We first consider the consequences of varying Ai while 
holding the other parameters fixed. For sufficiently small 



At, the quasienergy spectrum faithfully reproduces all 
the structure of the energy spectrum of the parent model. 
This is unsurprising, because if is greater than the 
bandwidth of the parent model's spectrum, direct multi- 
photon processes will not take place. If we now tune ujh 
so that it is less than this bandwidth, the quasienergy 
spectrum begins to deviate from the parent model's spec- 
trum at its edges. This effect can be seen, for instance, by 
examining the trace for At = 1 in panels (a) or (b). For 
even higher values of At (i.e. lower values of ojh), multi- 
photon processes strongly mix the states of the parent 
model, resulting in a flat quasienergy spectrum. 

The effect of multi-photon processes can also be en- 
hanced by broadening the parent model's spectrum, 
which can be achieved by raising g or u. In panel (c) of 
Figure 10 for instance, multi-photon processes have sig- 
nificantly flattened the spectrum for At = 1, and devia- 
tions from the parent model are even visible for At = 0.5. 
Since we always use At = 1 in our real-time dynamics 
simulations, it is perhaps fortunate that g = 0.9 is well 
within the proposed ergodic phase for u = 0.16 and that, 
near the critical point (i.e. in panel (b)), the quasienergy 
spectrum for At = 1 still retains much of the structure 
of the parent model's spectrum. 

However, there is one more caveat to keep in mind: 
the energy content of the system also grows with L. At 
flxed g, M, and At, the properties of the parent and 
time-dependent models deviate from one another as the 
system size grows. If we truly want to faithfully repro- 
duce the dynamics of the parent model with the modified 
dynamics, it may be necessary to scale At down as we 
raise L. However, recall that our goal is simply to find 
MBL in a model qualitatively similar to the parent model 
(1). Even with this more modest goal in mind, there is 
still the danger that, on sufficiently large lattices, multi- 
photon processes might couple a very large number of 
localized states and thereby destroy the many-body lo- 
calized phase of the parent model. Our numerical obser- 
vations indicate that this does not happen for the sys- 
tem sizes that we can simulate. We can keep At fixed at 
unity for L < 20 without issues, accepting the possibility 
that the sequence of models that we would in principle 
simulate on still larger lattices may require progressively 
smaller values of At. 



2. Level Statistics of the Many-Body Parent Model 

Localization transitions are often characterized by 
transitions in the level statistics of the energy spec- 
trum. Two of us previously looked at the level statis- 
tics of the disordered problem and identified a crossover 
from Poisson statistics in the many-body localized phase 
to Wigner-Dyson statistics in the many-body ergodic 
phase^"'^. The intuition that underlies this crossover is 
the following: in a localized phase, particle configurations 
that have similar potential energy are too far apart in 
configuration space to be efficiently mixed by the kinetic 



15 



0.4 
0.35 





















-s-0 
-a-0.04 
^0.16 
—0.32 








—0.64 













0.5 1 1.5 

9 



FIG. 11: The mean of the ratio between adjacent gaps in 
the spectrum, defined in (A4). This data was obtained by 
exact diagonahzation of the parent model (1) for L — 12 
systems. All data points have been averaged over 50 samples, 
and the legend refers to different values of the interaction 
strength u. The mean value of (r„) shows a crossover from 
Poisson statistics (indicated by the bottom reference line) to 
Wigner-Dyson statistics (indicated by the top reference line), 
for the largest values of u. Representative error bars have 
been included in the plot; the absent error bars have roughly 
the same size. 

energy term in the Hamiltonian. Therefore, level repul- 
sion is strongly suppressed, and Poisson statistics hold. 
Conversely, in an ergodic phase, there is strong level re- 
pulsion which lifts degeneracies, leading to Wigner-Dyson 
(i.e. random matrix) statistics. 

Along the lines of the aforementioned study of the dis- 
ordered problem, we focus on the gaps between successive 
eigenstates of the spectrum of the many-body "parent" 
model (1): 

= En+1 — En (A3) 



and a dimensionless parameter that captures the corre- 
lations between successive gaps in the spectrum: 

. --g4±4 (A4) 
max(c>„,d„+i) 

For a Poisson spectrum, the r„ are distributed as (^^^^^ 
with mean 2 ln(2) — 1 « 0.386; meanwhile, when ran- 
dom matrix statistics hold, the mean value of r has been 
numerically determined to be approximately 0.5295 ± 
0.0006". 

In Figure 11, we present exact diagonahzation results 
for L = 12 lattices at half-filling with potential wavenum- 
ber k — ^ and the boundary conditions described in Sec- 
tion II. C above. We show data for the same parameter 
range examined in the body of this paper and average 
over 50 samples for each value of g and u. For the largest 
value of u, the mean value of r„ interpolates between the 
expected values as g is raised, consistent with the exis- 
tence of a localization transition. We have also checked 
that the distributions of r„ have the expected forms in 
the small and large g limits in this regime. For smaller 
values of u, we can speculate that (r„) grows with L at 
large g and approaches the expected value for very large 
L. To argue for a MBL transition on the basis of exact 
diagonahzation, we would need to study this sharpening 
of the crossover as L is raised. This would indeed be an 
interesting avenue for future work. For our present pur- 
poses however, we only want to check consistency with 
our real-time dynamics data, as we have done in Figure 
11. 



1 P. Anderson, Physical Review 109, 1492 (1958). 

^ E. Abrahams, P. Anderson, D. Licciardello, and T. Ra- 

makrishnan. Physical Review Letters 42, 673 (1979). 
^ L. Fleishman and P. Anderson, Physical Review B 21, 2366 

(1980). 

* D. Thouless and S. Kirkpatrick, Journal of Physics C: Solid 
State Physics 14, 235 (1981). 

^ D. Basko, I. Aleiner, and B. Altshuler, Annals of physics 
321, 1126 (2006). 

® D. Basko, I. Aleiner, and B. Altshuler, Problems of Con- 
densed Matter Physics 1, 50 (2007). 

S. Sachdev, Quantum phase transitions (Wiley Online Li- 
brary, 2007). 

® V. Dobrosavljevic, Arxiv preprint arXiv:1112.6166 (2011). 
® J. Deutsch, Physical Review A 43, 2046 (1991). 
^° M. Srednicki, Physical Review E 50, 888 (1994). 

V. Oganesyan and D. Huse, Physical Review B 75, 155111 
(2007). 

A. Pal and D. Huse, Physical Review B 82, 174411 (2010). 
M. Znidaric, T. c. v. Prosen, and P. Prelovsek, Phys. Rev. 
B 77, 064426 (2008), URL http: //link . aps . org/doi/10 . 
1103/PhysRevB . 77 . 064426. 

A. Karahalios, A. Metavitsiadis, X. Zotos, A. Gorczyca, 



and P. Prelovsek, Phys. Rev. B 79, 024425 (2009), 
URL http : //link . aps . org/doi/10 . 1103/PhysRevB . 79 . 
024425. 

V. Oganesyan, A. Pal, and D. A. Huse, Physical Review B 
80, 115104 (2009). 
^'^ C. Monthus and T. Garel, Physical Review B 81, 134202 
(2010). 

^'^ T. Berkelbach and D. Reichman, Physical Review B 81, 
224429 (2010). 

E. Canovi, D. Rossini, R. Fazio, G. Santoro, and A. Silva, 
Physical Review B 83, 094431 (2011). 
J. Bardarson, F. Pollmann, and J. Moore, Arxiv preprint 
arXiv:1202.5532 (2012). 
^° R. Vosk and E. Ahman, Arxiv preprint arXiv: 1205.0026 
(2012). 

A. De Luca and A. Scardicchio, Arxiv preprint 
arXiv: 1206.2342 (2012). 
^2 S. Aubry and G. Andre, Arm. Israel Phys. Soc 3, 1 (1980). 
P. Harper, Proceedings of the Physical Society. Section A 
68, 874 (1955). 

D. Thouless and Q. Niu, Journal of Physics A: Mathemat- 
ical and General 16, 1911 (1983). 

J. Chaves and I. Satija, Physical Review B 55, 14076 



16 



(1997). 

^® S. N. Evangclou and D. E. Katsanos, Physical Review B 
56, 12797 (1997). 

^'^ A. G. Abanov, J. C. Talstra, and P. B. Wiegmann, Physical 
Review Letters 81, 2112 (1998). 

^® A. Eilmes, U. Grimm, R. Roemer, and M. Schreiber, The 
European Physical Journal B-Condensed Matter and Com- 
plex Systems 8, 547 (1999). 

L. Fallani, J. Lye, V. Guarrera, C. Fort, and M. Inguscio, 
Physical Review Letters 98, 130404 (2007). 
E. Lucioni, B. Deissler, L. Tanzi, G. Roati, M. Zaccanti, 
M. Modugno, M. Larcher, F. Dalfovo, M. Inguscio, and 
G. Modugno, Physical Review Letters 106, 230403 (2011). 
Y. Lahini, R. Pugatch, F. Pozzi, M. Sorel, R. Morandotti, 
N. Davidson, and Y. Silberberg, Physical Review Letters 
103, 13901 (2009). 

32 M. Albert and P. Leboeuf, Physical Review A 81, 013614 

(2010) . 

33 J. Cestari, A. Foerster, M. Gusmao, and M. Continentino, 
Physical Review A 84, 055601 (2011). 

3" N. Nessi and A. lucci. Physical Review A 84, 063614 

(2011) . 

35 ]y/j Tezuka and A. M. Garci'a-Garcfa, Physical Review A 

85, 031602 (2012). 
3'5 K. He, 1. Satija, C. Clark, A. Rey, and M. Rigol, Arxiv 

preprint arXiv:1201.2740 (2012). 
3^^ P. Ribeiro, M. Haque, and A. Lazarides, ArXiv preprint 

arXiv:1211.6012 (2012). 
3* C. Gramsch and M. Rigol, Arxiv preprint arXiv:1206.3570 

(2012) . 

3® C. Kittel and H. Kroemer, Thermal physics (WH Freeman 
& Co, 1980). 

*° M. Maricq, Physical Review B 25, 6622 (1982). 

T. Kitagawa, T. Oka, and E. Demler, Arxiv preprint 
arXiv:1201.0521 (2012). 

D. Martinez and R. Molina, The European Physical Jour- 
nal B-Condensed Matter and Complex Systems 52, 281 
(2006). 

*3 L. D'Alessio and A. Polkovnikov, ArXiv preprint 
arXiv:1210.2791 (2012). 

A. Kolovsky and G. Mantica, Arxiv preprint 
arXiv:1205.5078 (2012). 

Y. Lahini, Y. Bromberg, D. Christodoulides, and Y. Sil- 
berberg, Physical Review Letters 105, 163905 (2010). 
D. Basko, 1. Aleirier, and B. Altshuler, Physical Review B 
76, 052203 (2007). 
^"^ Both the many-body ergodic and localized phases differ 
qualitatively from their counterparts in the non-interacting 
AA model. The non-interacting extended phase is not er- 
godic, indicating that interactions are necessary for ther- 
malization. Meanwhile, the many-body localized phase 
is expected to exhibit logarithmic growth of the bipar- 
tite entanglement entropy to an extensive value, albeit 
with subthermal entropy density. Such behavior is in fact 
consistent with the recent observations in the disordered 
problem^3.i9 This growth is absent in the AA localized 
phase without interactions. Despite this difference, the in- 



teracting and non-interacting localized phases are similar 
in their inability to thermalize the particle density. 
There is an exception to this statement: multi-photon pro- 
cesses do seem to play an important role deep in the er- 
godic phase, where the energy content of the system is 
especially high. See Appendix A and the discussion of the 
time-dependence of the autocorrelator x in Section III. A 
for more details. 

'^^ To appropriately realize open boundary conditions, we 
should also turn off interactions over the boundary. When 
exploring different options for the boundary conditions, we 
varied J over the boundary and neglected to vary V. This is 
unfortunate in that it makes the model somewhat stranger. 
However, our boundary conditions are chosen for conve- 
nience anyway, and the numerics suggest that the choice of 
boundary conditions does not impact the essential physics 
discussed in this paper. 

^° The expectation value of the total energy of the many 
body system in a random state is of order ~ VL. Sup- 
pose this energy is conserved by the dynamics. We can 
write E/V~L = a;o -f hAo cos 6*0 ~ Xoc + hAoc cos 9aa- Here, 
the subscripts and oo refer to the initial and late-time 
states, xo and bounded random numbers capturing 

the expectation value of interactions (and hopping at late- 
times) , h is the non-random amplitude of the quasiperiodic 
potential, and Ao and Aac are positive bounded ampli- 
tudes of the Fourier components at the wavevector k of the 
quasiperiodic potential. This ansatz implies a finite correla- 
tion between the random phases 6q and 6aa ■ Therefore, one 
of the Fourier modes of x remains correlated as L — > oo, 
and we expect x ~ ;^ in the ergodic phase. Note that this 
argument truly applies only to the energy-conserving par- 
ent model. In fact, in our rmmcrics, there is only partial 
energy conservation, and energy non-conserving events be- 
come more prevalent as u, g, or L is raised. This means 
that X will generically decay faster than at large L in 
the ergodic phase. 

This statement should be interpreted with some care. 
Quantum entanglement entropy measures, such as the 
Rcnyi entropy that we define in equation (14), carry in- 
formation about the off-diagonal elements in the reduced 
density matrix. These terms have no classical analogue and 
would not be considered in a thermodynamic calculation. 
This difference can result in discrepancies in the sublead- 
ing behavior. For instance, consider our calculation of the 
bipartite Rcnyi entropy of the model state |$) in Section 
IV. A: the quantum Renyi entropy is one bit lower than 
the Renyi entropy calculated by classical counting of con- 
figurations. A more precise analogue of the classical en- 
tropy would thus be a "diagonal" entropy in which all 
off-diagonal elements of the reduced density matrix were 
neglected. 

*2 Only the first term on the right-hand side of equation (22) 
would appear in a "classical counting" derivation of the 
thermodynamic entropy. The other two terms account for 
off-diagonal elements in the reduced density matrix (20). 
Please see footnote 51 for more details. 



